1 Introducción

En esta práctica utilizaremos los registros de la especie Leopardus wiedii, la paquetería ENMEval 2.0 (Kass et. al 2025) y el algoritmo MaxEnt (Phillips et al. 2006) para explorar diferentes métodos de partición de datos en la modelación de distribución de especies. El objetivo de esta práctica es comparar tipos de partición de datos, evaluando cómo cada enfoque maneja la autocorrelación espacial entre los datos de entrenamiento y validación, así como su efecto en las métricas.

2 Calibración

2.1 Carga de datos biológicos y ambientales

# Leer datos biológicos
datos <- read.csv("/Users/jimenaburgos/Downloads/practica_8p1/train_lw.csv")
occs <- datos[,c("long","lat")]

# Leer los datos ambientales
setwd("/Users/jimenaburgos/Downloads/CursoNichos_2026-1/Practica2_Evaluation/envs")
envs.files <- list.files(pattern='.tif', full.names=TRUE)

# Hacemos el raster stack 
envs <- rast(envs.files)

# Graficamos para ver los registros
plot(envs[[3]])
points(occs, pch=19)

2.2 Carga de capa de sesgo

### Cargamos una capa de sesgo para dirigir los puntos de background
biasfile = rast("/Users/jimenaburgos/Downloads/CursoNichos_2026-1/Practica2_Evaluation/biaslayer_HREFBiv.tif")
biaslayer <- biasfile
# Es necesario que no existan NAs en el raster
# Por lo que asignamos un valor de 0 a los NAs
values(biasfile)[is.na(values(biasfile))] <- 0
n.bg <- 1000 # este sera la cantidad de puntos de background 
bg <- xyFromCell(biasfile, 
                 sample(ncell(biasfile),
                        n.bg, prob=values(biasfile))) %>% as.data.frame()
colnames(bg) <- colnames(occs)
head(bg)
##        long      lat
## 1 -87.02083 20.43750
## 2 -90.43750 19.89583
## 3 -87.81250 18.97917
## 4 -87.89583 20.81250
## 5 -91.06250 18.18750
## 6 -87.39583 21.22917
# Graficamos los puntos de background sobre la capa de sesgo
plot(biaslayer, main = "Capa de sesgo")
points(bg, pch=3 , cex = 0.4, col="white")

2.3 Función pROC

# Definimos una función personalizada para implementar la pROC
proc <- function(vars) {
  proc <- kuenm::kuenm_proc(vars$occs.val.pred, c(vars$bg.train.pred, vars$bg.val.pred))
  out <- data.frame(proc_auc_ratio = proc$pROC_summary[1], 
                    proc_pval = proc$pROC_summary[2], row.names = NULL)
  return(out)
}

3 Particiones

3.1 Particiones aleatorias

Los dos métodos siguientes no tienen en cuenta la autocorrelación espacial entre los registros de validación y entrenamiento. Por lo que, es importante considerar que la partición aleatoria puede llevar a la sobreestimación del modelo si los registros de presencia de calibración y evaluación están próximas entre sí, debido a que las localidades utilizadas para evaluar el modelo no son independientes de las utilizadas para calibrarlo.

3.1.1 Jackknife (leave-one-out)

Principalmente, al trabajar con conjuntos de datos relativamente pequeños (e.g. 25 registros de presencia), los usuarios pueden optar la validación cruzada k-fold, en el que el número de bins (k) o grupos es igual al número de registros de presencia (n) (Pearson et al., 2007; Shcheglovitova y Anderson, 2013). Esto se conoce como partición jackknife o “dejar uno fuera” (Hastie et al., 2009).

jack <- get.jackknife(occs, bg)
evalplot.grps(pts = occs, pts.grp = jack$occs.grp, envs = envs)

# Corremos el modelo
e.mx_jack <- ENMevaluate(occs = occs, #Registros de presencia
                          envs = envs, #Capas ambientales
                          bg = bg, #Puntos de background
                          algorithm = 'maxnet', #Algoritmo
                          partitions = 'jackknife', #Método de partición
                          tune.args = list(fc = c("L","LQ","LQP"), rm = 1),
                          user.eval = proc) #Métrica adicional (pROC)
##   |                                                                              |                                                                      |   0%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================================| 100%

3.1.2 Random k-fold

El método de «k-fold aleatorio» divide aleatoriamente los registros de presencia en un número de (k) bins o grupos especificado por el usuario (Hastie et al., 2009). Este método es ideal utilizarlo con conjuntos de registros de presencia grandes.

rand <- get.randomkfold(occs, bg, k = 5) #En este caso elegimos 5
evalplot.grps(pts = occs, pts.grp = rand$occs.grp, envs = envs)

# Corremos el modelo
e.mx_kfold <- ENMevaluate(occs = occs, 
                        envs = envs,
                        bg = bg, 
                        algorithm = 'maxnet', 
                        partitions = 'randomkfold', 
                        tune.args = list(fc = c("L","LQ","LQP"), rm = 1),
                        user.eval = proc)
##   |                                                                              |                                                                      |   0%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================================| 100%

3.2 Particiones geográficamente estructuradas

Los métodos de partición son variaciones de las particiones geográficamente estructuradas (Radosavljevic y Anderson, 2014). Básicamente, estos métodos dividen los registros de presencia y background en grupos de evaluación basados en reglas espaciales. El objetivo es reducir la autocorrelación espacial entre los registros utilizados en la validación y el entrenamiento, evitando inflar el rendimiento del modelo, al menos para los conjuntos de datos resultantes de un muestreo sesgado (Veloz 2009, Wenger y Olden 2012, Roberts et al. 2017).

3.2.1 Bloques

# El método de "bloque" particiona los datos según la mediana de la latitud 
# y longitud que dividen los registros de presencia en cuatro grupos espaciales de
# igual número o lo más próximos posible.

block <- get.block(occs, bg, orientation = "lat_lon")

# El objeto resultante es una lista de dos vectores que proporcionan la
# designación del bloque o bin para cada registro de presencia y background
table(block$occs.grp)
## 
##  1  2  3  4 
## 10 10 10  9
# La primera dirección divide los registros de presencia en dos grupos, y la 
# segunda divide cada uno de ellos en dos grupos adicionales, resultando en cuatro grupos.
# Tanto los registros de presencia como los de background se asignan a cada uno de 
# los cuatro bloques o bins, según su posición con respecto a estas líneas 

# bloques de registros de presencia
evalplot.grps(pts = occs, pts.grp = block$occs.grp, envs = envs) + 
  ggplot2::ggtitle("Spatial block partitions: occurrences")

# bloques de background
evalplot.grps(pts = bg, pts.grp = block$bg.grp, envs = envs) + 
  ggplot2::ggtitle("Spatial block partitions: background")

# Corremos el modelo
e.mx_block <- ENMevaluate(occs = occs, 
                     envs = envs,
                     bg = bg, 
                     algorithm = 'maxnet', 
                     partitions = 'block', 
                     tune.args = list(fc = c("L","LQ","LQP"), rm = 1),
                     user.eval = proc)
##   |                                                                              |                                                                      |   0%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================================| 100%

3.2.2 Checkerboard básico

Estos generan cuadrículas como un tablero de ajedrez a lo largo de la extensión del área de calibración y dividen los registros de presencia en grupos según su ubicación en el tablero. A diferencia del método de bloques, ambos métodos de tablero de ajedrez subdividen el área de calibración equitativamente, aunque no garantizan el mismo número de localidades de ocurrencia en cada bin o grupo.

El método de tablero de ajedrez básico divide los puntos en k = 2 grupos utilizando un patrón de tablero de ajedrez simple.

cb1 <- get.checkerboard(occs, envs, bg,
                         aggregation.factor=30) #puedes cambiar
                          #este parámetro para crear "cajas" más grandes.
# visualizamos presencias
evalplot.grps(pts = occs, pts.grp = cb1$occs.grp, envs = envs)

# visualizamos background
evalplot.grps(pts = bg, pts.grp = cb1$bg.grp, envs = envs)

# Corremos el modelo
e.mx_check1 <- ENMevaluate(occs = occs, 
                         envs = envs,
                         bg = bg, 
                         algorithm = 'maxnet', 
                         tune.args = list(fc = c("L","LQ","LQP"), rm = 1),
                         partitions = "user",
                         user.grp = cb1,
                         user.eval = proc)
##   |                                                                              |                                                                      |   0%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================================| 100%

3.2.3 Checkerboard jerárquico

El método de tablero de Checkerboard jerárquico genera k = 4 grupos espaciales mediante un enfoque de dos escalas anidadas. Este método aplica dos factores de agregación secuenciales: primero crea una grilla gruesa que divide el área de estudio, y luego subdivide cada celda de esta grilla en celdas más finas, generando un patrón espacial jerárquico.

A diferencia del checkerboard básico que produce solo 2 grupos, este método ofrece una partición más granular con 4 grupos, manteniendo la separación geográfica pero proporcionando mayor flexibilidad para conjuntos de datos grandes o áreas de estudio complejas. Los registros de presencia y background se asignan a cada grupo según su ubicación en el patrón jerárquico resultante.

cb2 <- get.checkerboard(occs, envs, bg, aggregation.factor=c(10,10))
# visualizamos presencias
evalplot.grps(pts = occs, pts.grp = cb2$occs.grp, envs = envs)

# visualizamos background
evalplot.grps(pts = bg, pts.grp = cb2$bg.grp, envs = envs)

# Corremos el modelo
e.mx_check2 <- ENMevaluate(occs = occs, 
                          envs = envs,
                          bg = bg, 
                          algorithm = 'maxnet', 
                          tune.args = list(fc = c("L","LQ","LQP"), rm = 1),
                          partitions = "user",
                          user.grp = cb2,
                          user.eval = proc)
##   |                                                                              |                                                                      |   0%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================================| 100%

4 Métricas

4.1 Resultados de evaluación

# Podemos ver los resultados de cualquiera de nuestros modelos
# vemos el de jackknife
e.mx_check1@results
##    fc rm   tune.args auc.train cbi.train auc.diff.avg auc.diff.sd auc.val.avg
## 1   L  1   fc.L_rm.1 0.6900187     0.771   0.06782160  0.07391116   0.6574025
## 2  LQ  1  fc.LQ_rm.1 0.6902389     0.895   0.06133692  0.06881793   0.6577346
## 3 LQP  1 fc.LQP_rm.1 0.6990200     0.901   0.06518888  0.08016369   0.6558617
##   auc.val.sd cbi.val.avg cbi.val.sd or.10p.avg or.10p.sd or.mtp.avg  or.mtp.sd
## 1 0.05742105      0.4685 0.01343503  0.2045455 0.2892710 0.11363636 0.16070609
## 2 0.05403649      0.4565 0.05161880  0.1818182 0.2571297 0.09090909 0.12856487
## 3 0.06439568      0.3720 0.20364675  0.2045455 0.2892710 0.06818182 0.09642365
##   proc_auc_ratio.avg proc_auc_ratio.sd proc_pval.avg proc_pval.sd     AICc
## 1           1.167737        0.07570549             0            0 897.1241
## 2           1.177875        0.05026016             0            0 900.6929
## 3           1.159173        0.08067267             0            0 902.5415
##   delta.AICc      w.AIC ncoef
## 1   0.000000 0.81002561     4
## 2   3.568711 0.13600752     5
## 3   5.417391 0.05396688     6

4.2 Compilación de resultados

# Asignamos los resultados de la evaluacion en un dataframe
res_block <- eval.results(e.mx_block)
res_check1 <- eval.results(e.mx_check1)
res_check2 <- eval.results(e.mx_check2)
res_jack <- eval.results(e.mx_jack)
res_kfold <- eval.results(e.mx_kfold)

4.3 Comparación y visualización

# Agregamos una columna para distinguir los tipos de partición
res_block["Tipo_particion"] <- "Bloque"
res_check1["Tipo_particion"] <- "CheckerboardBasico"
res_check2["Tipo_particion"] <- "CheckerboardJerarquico"
res_jack["Tipo_particion"] <- "Jackknife"
res_kfold["Tipo_particion"] <- "Random-k-fold"

# Unimos los dataframes
data_results <- rbind(res_block, res_check1, res_check2, res_jack, res_kfold)

data_long <- data_results %>%
  # Seleccionar las columnas necesarias
  select(Tipo_particion, or.10p.avg, auc.val.avg ) %>%
  # Convertir a formato largo directamente (sin agrupar ni promediar)
  pivot_longer(
    cols = c(or.10p.avg, auc.val.avg),
    names_to = "Metrica",
    values_to = "Valor"
  ) %>%
  # Renombrar las métricas para que se vean mejor en la gráfica
  mutate(
    Metrica = case_when(
      Metrica == "or.10p.avg" ~ "OR", 
      Metrica == "auc.val.avg" ~ "AUC",
      TRUE ~ Metrica
    )
  )
# Visualizamos las métricas de los distintos tipos de partición
ggplot(data_long, aes(x = Tipo_particion, y = Valor, color = Tipo_particion, fill = Tipo_particion)) +
  geom_boxplot(alpha = 0.7) +
  facet_wrap(~Metrica, scales = "free_y") +
  theme_minimal() +
  labs(
       x = "Tipo de partición",
       y = "Valor") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none")

5 Visualización de modelos

5.1 Selección de modelos óptimos

# Utilizamos un protocolo secuencial para elegir un modelo de cada tipo de partición:
# la menor tasa de omisión y el máximo AUC
opt.seq_block <- res_block %>% 
  filter(or.10p.avg == min(or.10p.avg)) %>% 
  filter(auc.val.avg == max(auc.val.avg))

opt.seq_check1 <- res_check1 %>% 
  filter(or.10p.avg == min(or.10p.avg)) %>% 
  filter(auc.val.avg== max(auc.val.avg))

opt.seq_check2 <- res_check2 %>% 
  filter(or.10p.avg == min(or.10p.avg)) %>% 
  filter(auc.val.avg == max(auc.val.avg))

opt.seq_jack <- res_jack %>% 
  filter(or.10p.avg == min(or.10p.avg)) %>% 
  filter(auc.val.avg == max(auc.val.avg))

opt.seq_kfold <- res_kfold %>% 
  filter(or.10p.avg == min(or.10p.avg)) %>% 
  filter(auc.val.avg == max(auc.val.avg))

5.2 Extracción de modelos

# Ahora seleccionamos solo ese modelo de nuestro conjunto de modelos
mod.seq_block <- eval.models(e.mx_block)[[opt.seq_block$tune.args]]
mod.seq_check1 <- eval.models(e.mx_check1)[[opt.seq_check1$tune.args]]
mod.seq_check2 <- eval.models(e.mx_check2)[[opt.seq_check2$tune.args]]
mod.seq_jack <- eval.models(e.mx_jack)[[opt.seq_jack$tune.args]]
mod.seq_kfold <- eval.models(e.mx_kfold)[[opt.seq_kfold$tune.args]]

5.3 Predicciones y mapas

# Finalmente seleccionamos una predicción
pred.seq_block <- eval.predictions(e.mx_block)[[as.character(opt.seq_block$tune.args)]]
pred.seq_check1 <- eval.predictions(e.mx_check1)[[as.character(opt.seq_check1$tune.args)]]
pred.seq_check2 <- eval.predictions(e.mx_check2)[[as.character(opt.seq_check2$tune.args)]]
pred.seq_jack <- eval.predictions(e.mx_jack)[[as.character(opt.seq_jack$tune.args)]]
pred.seq_kfold <- eval.predictions(e.mx_kfold)[[as.character(opt.seq_kfold$tune.args)]]
# Graficamos los modelos seleccionados de cada tipo de partición
par(mfrow=c(1, 3))
plot(pred.seq_block, main="Block")
plot(pred.seq_kfold, main="Random-k-fold")
plot(pred.seq_check1, main="Checkerboard básico")

par(mfrow=c(1, 2))
plot(pred.seq_check2, main="Checkerboard jerárquico")
plot(pred.seq_jack, main="Jackknife")